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A Lagrangian relativistic approach to the non-linear dynamics of cosmological pertur- 
bations of an irrotational collisionless fluid is considered. Solutions are given at second 
order in perturbation theory for the relevant fluid and geometric quantities and com- 
pared with the corresponding ones in the Newtonian approximation. Specifically, we 
compute the density, the volume expansion scalar, the shear, the "electric" part, or 
J> ■ tide, and the "magnetic" part of the Weyl tensor. The evolution of the shear and the 

, tide beyond the linear regime strongly depends on the ratio of the characteristic size 

of the perturbation to the cosmological horizon distance. For perturbations on sub- 
, horizon scales the usual Newtonian approximation applies, at least at the considered 

perturbative order; on super-horizon scales, instead, a new picture emerges, which we 
t^J" , call "silent universe", as each fluid element evolves independently of the environment, 

0^ ' being unable to exchange signals with the surrounding matter through either sound 

^ <— | waves or gravitational radiation. For perturbations inside the Hubble radius particu- 

lar attention is paid in singling out non-local effects during the non-linear evolution 
of fluid elements. These non-local effects are shown to be carried by a traceless and 
■ divergenceless tensor, contained in the magnetic part of the Weyl tensor, which is 

, dynamically generated as soon as the system evolves away from the linear regime. 
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1 INTRODUCTION 

The study of the dynamics of self-gravitating fluids is of extreme importance in the study of the formation of structures 
in the universe for from primordial perturbations via the gravitational instability. The solutions of the relevant dynamical 
equations are obtained by different techniques depending on the specific application. A Newtonian approach is usually applied 
for perturbations on scales much smaller than the horizon distance, but much larger than the Schwarzschild radii of collapsing 



objects (Peebles 1980). Analytical methods are then used to follow the evolution of small perturbations and useful approximate 
solutions are also available to describe the system up to the mildly non-linear regime, such as the celebrated Zel'dovich 
approximation, while the general non-linear dynamics can only be followed numerically, e.g. by N-body codes. A relativistic 
treatment in cosmology is instead applied to study linear perturbations on large scales and/or at early times, and to the 
analysis of those few symmetric solutions of Einstein's equations which are relevant in a cosmological context (e.g. the 
spherically symmetric Tolman-Bondi metric). 

Having a unified treatment of the problem, able to describe both the large-scale properties of the universe and the non- 
linear collapse of structures by gravitational instability on small scales would be highly desirable. Building up such a theory 
is however hampered by the absence of symmetries in the system, which implies that all the matter and geometric degrees of 
freedom are possibly switched on during the non-linear evolution of self-gravitating matter. 

A possible picture of the problem is represented by the fluid-flow approach, in which the dynamics is followed directly 
in terms of observable fluid variables (such as density, shear, vorticity, etc.) and other tensor quantities which describe the 
space-time curvature (instead of the metric tensor). A complete description of the method can be found in the classical 
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review by Ellis (1971). Following this line of thought, Matarrese, Pantano & Saez (1993) have recently proposed a Lagrangian, 
general relativistic approach to the fully non-linear dynamics of a collisionless fluid. The approach is based on two relevant 
simplifying assumptions: vanishing vorticity (i.e. irrotational flow) and negligible signal exchange by the gravitational field, as 
it is obtained by dropping the terms containing the so-called magnetic part H a b of the Weyl tensor in the non-linear evolution 
of the tide (or electric part of the Weyl tensor, E a b)- While the former assumption is a widely accepted one in the standard 
gravitational instability scenario of structure formation, the latter is certainly more problematic. The meaning of H a b in linear 
theory is well known: it only contains vector and tensor modes (e.g. Bruni, Dunsby & Ellis 1992); if the vorticity vanishes, 
as in our case, no physical vector modes are present and H a b at most contains gravitational waves. Beyond linear theory, 
however, the physical content of H a b is less clear. One can reasonably assume that in the absence of H a b no gravitational 
waves occur: this is evident in a Lagrangian approach, where, the absence of pressure gradients and the vanishing of H a b 
imply that no spatial gradients appear in the evolution equations (apart from those possibly contained in convective time 
derivatives, which can be dropped by going to the comoving matter frame). It is then clear that no actual wave propagation 
can occur when spatial derivatives are absent in the fluid and gravitational evolution equations. The method proposed by 
Matarrese, Pantano & Saez (1993) takes advantage of the absence of spatial gradients to construct a Lagrangian algorithm in 
which the dynamics of fluid elements is completely determined by local initial conditions; each fluid element can be therefore 
independently followed up to the time when it enters a multi-stream region. In this sense the method is quite similar to the 
well-known Zel'dovich approximation in the Newtonian context, to which it actually reduces to first order in a perturbation 
expansion. 

While neglecting the fluid pressure (i.e. taking p — 0) and considering zero initial vorticity Lu a b are reasonable assumptions 
in a cosmological framework, the condition H a b = has to be taken with more caution. It has been shown (Barnes & 
Rowlirgson 198E) that the only solutions of Einstein equations, with p = uo a b = H a b = are either of Petrov type I, or 
conformally flat, or homogeneous and anisotropic of Bianchi type I, or with two degenerate shear eigenvalues and described 
by a Szekeres line-element. All of these cases require some restrictions on the initial data, which means that the exact 
conditions above are not suitable for studying cosmological structure formation with enough generality; requiring p = ui a b = 
and Hab ~ appears more feasible. In a recent work Matarrese, Pantano & Saez (1994) have clarified the dynamical role 
of the magnetic tensor, by calculating the behaviour of fluctuations around Robertson- Walker (RW) at second order in 
perturbation theory: whatever initially scalar perturbations are given, H a b vanishes at first order, but not beyond. The 
important consequence of this result is that the non-vanishing H a b allows for the influence of the surrounding matter on the 
evolution of fluid elements, through the presence of spatial gradients in the evolution equations. The magnetic tensor appears 
therefore responsible for the environmental influence of the matter on the non-linear evolution of fluid elements. As we shall 
see below, although the signal carrying such a non-local information travels at finite speed, for perturbations on scales much 
smaller than the horizon distance it appears as an instantaneous Newtonian effect. 

In this paper we give a detailed derivation of the perturbative solutions at second order of the fluid and gravitational field 
equations, expressed in Lagrangian form. We pay particular attention to the occurrence of non-local effects in the evolution of 
fluid elements away from the linear regime. This work complements and extends the results previously reported in (Matarrese, 
Pantano & Saez 1994). 

The plan of the paper is as follows. In Section 2 we review the general fluid-flow method used to study the relativistic 
dynamics of a self-gravitating pressureless fluid, which we then specialize to the case of irrotational motions represented in a 
comoving frame, i.e. in Lagrangian form; to this aim we introduce suitably rescaled variables which describe deviations out of 
a homogeneous and isotropic universe expansion. In Section 3 we present our second order expansion and discuss its properties 
on scales larger and smaller than the horizon distance. Some cosmological consequences of our work are discussed in Section 
4. Two appendices are added to help the comparison with the Newtonian approximation and to provide some technical details 
on the derivation of the second order expressions. 

We use the signature +2; latin indices refer to space-time coordinates, (1,2,3,4), greek indices to spatial ones, (1,2,3); 
semicolons are used for covariant derivatives, commas for partial ones. Units are such that the speed of light c is one. 



2 RELATIVISTIC DYNAMICS 

In this Section we shall introduce the equations which govern the dynamics of a collisionless perfect fluid in general relativity. 
The dynamical equations are written in terms of observable fluid quantities and other tensor quantities which directly describe 
the space-time curvature. A complete description and a full derivation of the equations presented here can be found in (Ellis 
1971 . 



© 0000 RAS, MNRAS 000, 000-000 



Relativistic Approach to Gravitational Instability in the Expanding Universe 3 



2.1 Relativistic hydrodynamical equations 

The relativistic dynamics of a collisionless (i.e. with vanishing pressure) self-gravitating perfect fluid is determined by Einstein's 
field equations and by the continuity equations for the matter stress-energy tensor T ab = gu a u b , where g is the energy density 
(or matter density, in this case) and u a the four-velocity of the fluid, normalized to u a u a = — 1. The local projection tensor 
into the rest frame of an observer moving with four-velocity u a is defined as h ab = g ab + u a u b , for which h ab u b = 0. By 
differentiating the velocity field one obtains the space-like tensor v ab = h a c h b d u c -d, v ab u b = 0, and the acceleration vector 
ii a = u a . b u b , which is also space-like, ii a u a — 0, as it follows from the u a normalization. An overdot denotes convective 
differentiation with respect to the proper time t of fluid elements (for a general n-index tensor, A aia2 ... an = A aia2 ... an - b u b ). 

Being flow-orthogonal, the tensor v ab can only have nine independent components: its skewsymmetric is the vorticity 
tensor uj ab = v\ a a (the symbol r i stands for skewsymmetrization, r \ for symmetrization) , which describes rigid rotations of 
fluid elements with respect to a locally inertial rest frame. Due to its antisymmetric character, u ab has only three independent 
components which correspond to the vorticity vector u> a = ^rj abcd u b ijj cd , with r\ abcd the four-index Levi-Civita tensor. Equiv- 
alently, uj ab = r) abc< iuj c u d . The symmetric part of v ab can be further split into its trace O = v a a , called the volume expansion 
scalar and its traceless part, a ab = v^ a b) — ^Qh ab , called the shear tensor. The volume expansion scalar, giving the local rate 
of isotropic expansion (or contraction), defines a length-scale £, through the equation O = 31/1. This length-scale is just the 
generalization of the scale-factor a(t) of homogeneous and isotropic RW models; in that particular case O = 3H, where H(t) 
is Hubble's constant. The shear, on the other hand, describes a pure straining in which a spherical fluid volume is distorted 
into an ellipsoid with axis lengths changing at rates determined by the three a a b eigenvalues, ai, 02 and 03 = — (<7i +02). The 
vanishing trace condition implies that this deformation leaves the fluid volume invariant, while, in the absence of vorticity, 
the principal axes of the shear keep their direction fixed during the evolution, in a local inertial rest frame. 

For a collisionless perfect fluid the energy-momentum local conservation equations reduce to 

ii a = 0. (1) 

q = -t?e, (2) 

Eq.(l) tells us that in absence of pressure gradients each fluid element moves along a geodesic, while Eq.(2) is equivalent to 
the continuity equation for the particle number density, since the matter and energy density are directly proportional in our 
case. 

The expansion scalar satisfies the Raychaudhuri equation 

9 = A - ie 2 + 2(co 2 - a 2 ) - 4nGg, (3) 
o 

where G is Newton's constant and we have introduced a possible cosmological constant A and the two scalars u> 2 = u] a uj a = 
^ui ab ui ab and a 2 = \a ab a a b- Note that in the homogeneous and isotropic RW case, u> ab = a ab — and the latter equation 
reduces to the familiar Friedmann equation, 3(H + H 2 ) — —AnGg + A. 
The vorticity vector evolves according to 
9 

uj — —-<Slo + a b cu , (4) 

while the shear is determined by the equation 
1 2 

o~ab = -o ac a c b - uj a (j b + -h ab (2a + oj ) - -Qa ab - E ab , (5) 

where E ac = C abc dU b u d is the electric part of the Weyl tensor C a bcd- The tensor E ab is also called tidal tensor, since it contains 
that component of the gravitational field which describes tidal interactions; it is symmetric, traceless and flow-orthogonal, 
E ab u b = 0. 

From the Weyl tensor one can define another quantity: its magnetic part H ac = \r\ ah 9h C g hcdU b u d which is also symmetric, 
traceless and flow-orthogonal. While the tidal field has a straightforward Newtonian analogue, which can be written in terms 
of derivatives of the gravitational potential, the magnetic tensor H ab has no straightforward interpretation in the Newtonian 
equations. The important point is that, while in Newtonian theory, in its standard formulation, the gravitational potential 
is determined through a constraint, namely Poisson's equation, in general relativity both E ab and H ab can be calculated by 
solving suitable evolution equations. 

From the Bianchi identities one can derive the following evolution equations for E ab and H ab 

Eab = —h c ( a h b)m r\ mrad u r H c a .d — h ab a cd E cd — QE ab + E c ( a u) b) c + 3E c ( a a b) c — 4iTGgcr ab , (6) 

Hab = h c (ahb) m ri mrsd u r E c s .d — h a b^ cd H cd - OH a b + ffc(a^) C + 3H a ( a & b ) C - (7) 
In a collisionless fluid, the velocity field satisfies the following constraint equations 

= 0. (8) 
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h a b (u 1 h %-a b %+^& b ) = 0, (9) 
Further constraints need to be satisfied by the electric and the magnetic tensor, namely 

H ad = h*h£{w£' c + cr (t f ' ;c )77 i ,) /6c M / (10) 

h t a E as i dh s d - rf bpq u b a p d H qd + 3i?' s a; s = ^^-h t b p ,b (11) 

^ a H aa ;dK - rf bpq u b a d E qd + 3E t s uj s = p<J (12) 

Finally, another useful quantity is the (infinitesimal) relative position vector x^ between two neighbouring particles 
defined in the rest frame of the observer moving with four-velocity u a (x_i_ a u a = 0). The vector components x® evolves 
according to the "generalized Hubble law" 

x± = ~6 xl + (uj a b + a\) x b . (13) 
2.2 Cosmological equations 

In a cosmological context one can safely assume zero initial vorticity for the fluid, i.e. u) a — 0; Kelvin's circulation theorem 
implies that this condition continues to hold until the fluid enters a multi-stream region. Under the assumption of vanishing 
magnetic tensor, Matarrese, Pantano & Saez (1993) proposed a method for integrating the evolution equations in local inertial 
frames, taking advantage of the geodesies motion of collisionless particles. However, when spatial gradients explicitly occur, 
as it happens when H ab 7^ 0, it is more convenient to adopt a global coordinate system, such as one comoving with the fluid. 
Let us then rewrite the equations of Section 2.1 in comoving coordinates. Recalling that for a pressureless fluid comoving 
hypersurfaces (orthogonal to the energy flow) and synchronous ones (orthogonal to geodesies) coincide, one can write the 
line-element in the form 

ds 2 = — dt 2 + a 2 (t)h a pdq a dq 13 , (14) 

where a = At 2 / 3 , as for a flat, matter-dominated RW model with vanishing cosmological constant (our "background" solu- 
tion) and the rescaled spatial projection tensor h a f) generally depends upon space and time coordinates. For computational 



convenience we then introduce suitably rescaled quantities ( Matarrese, Pantano fc Saez 1994 ): a scaled density fluctuation 
A = (6-KGt 2 g — l)/o, a peculiar expansion scalar 1? = (3t/2a)(0 — 2/t), a traceless shear tensor s"^ = (3t/2a)a a g and a 
traceless tidal tensor e a g = (3t 2 /2a)E a g. These quantities can be grouped in two space-like tensors: the velocity gradient 
tensor 1?^ = s a p + |i5 a g#, related to the covariant derivatives of the peculiar velocity field and the peculiar gravitational field 
tensor A a p = e a p + ^AS^. We also rescale the magnetic tensor as H a p = (3t 2 /2a)H a /3 . The dynamical equations for the fluid 
and the gravitational field then read 



. a I 



+^h 0v (r)^ s H\, s +fj a - yS H\,s) , (16) 

ny = --n^ - 2m% - s a ^ s H\ + \n a ^ l3 + ^tc p ^ - ^^(^a% ;5 +fj^ s A\, s ) , (17) 

where the prime denotes partial differentiation with respect to the scale factor a and rj a ~' s is the Levi-Civita tensor relative to 
the conformal spatial metric h a /3, namely rj a/31 = h~ 1 ^ 2 e al}l , with e 123 = 1. The traces of Eqs.(15) and (16) yield, respectively, 

tf' = ~(tf + A)-tf a 7 ^ a) (18) 

A' = -&A- -0& + A), (19) 
a 

while the trace of Eq.(17) vanishes identically. As shown by Matarrese, Pantano & Saez (1994), Eq.(15) coincides with the 
equation obtained by taking the spatial gradients of the Newtonian Euler equation in Lagrangian form; Eq.(19) is instead 
identical to the Newtonian continuity equation, also in Lagrangian form (see Appendix Al). Since the traceless component of 
the evolution equation for A^ (i.e. the evolution equation for the tide) has no Newtonian analog, it follows that the magnetic 
tensor which enters through its spatial gradients in that equation, has no direct dynamical role in the Newtonian theory as it 
is usually formulated (i.e. in the standard form of Appendix Al). 

In our comoving gauge, the spatial metric tensor, which appears explicitly and through its spatial gradients in the 
covariant derivatives of Eqs.(16) and (17), evolves according to the equation 
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h^h'^ = #V (20) 
Using now our rescaled variables we can rewrite the constraints Eqs.(9)-(12) as 

tfcf;/3 = tf, Q , (21) 

Hi = ±h a ,(r jS ^; s +ff'*f', s ), (22) 

A£; = A, Q ~h ail hp v ff^r x H p ^ (23) 

H a % = ^hcjip^r^. (24) 

Note that the only Christoffel symbols involved in the covariant derivatives of Eqs.(16), (17) and (21) - (24) are the 
spatial ones, i.e. = T^ 7 , which can be trivially obtained from h a 0- 
The rescaled displacement vector = x" /a obeys the equation 

r' = *V- (25) 

The matrix connecting the Eulerian coordinates x a with the Lagrangian ones (f is the Jacobian 

J a g = dx a /dq = S a B +V% (26) 

where T> a j 3 is the symmetric deformation tensor. Taking then f a = dx a = J"*p£(Q\, where {J^ = dq 13 represents the initial (i.e. 
Lagrangian) infinitesimal displacement, one gets 

v y = 0- p + <f*fPp (27) 

which is formally solved by 



X>°s(a) = exp / da d%a) - 5% (28) 

Jag 

Once the Jacobian is known one gets the metric at the "time" a by performing the transformation 
h a p(a) = h J s(ao)f y a J p- (29) 

2.3 Initial conditions 

Initial perturbations can be given at the "time" ao in a RW bac kground by using the gauge-invar iant formalism (Bardeen 
1980). Here we consider only growing-mode scalar perturbations ( Vlatarrese, Pantano fc Saez 199sj ), 



A a p(a ) = -# Q a(ao) = ^o,"/3 , (30) 

where the scalar ipo, an arbitrary function of the spatial coordinates q a , corresponds to the initial peculiar gravitational 
potential, which is related to Bardeen's gauge-invariant $h by </?o = — (3/2A 3 )$h- These initial conditions correspond to the 
"seed" metric 

20 

h a p = S a g(l —A ip ) — 2a(p , a f3 , (31) 

and imply vanishing magnetic tensor Ti. a p at the linear level. In what follows, we will only consider the growing mode, 
proportional to atpo, a p, since the constant mode, proportional to A 3 ipo <C 1, can be practically neglected during the subsequent 
evolution. Also, in many cases it will be convenient to assume that the initial conditions are given at ao — > 0. 



It is important to realize that the constraints, Eqs.(21)-(24), are all satisfied at the linear level (Bruni, Dunsby & Ellis 



1992) by our choice of initial conditions. 



3 SECOND ORDER SOLUTIONS: BEYOND THE ZEL'DOVICH APPROXIMATION 

We shall now construct a second order Lagrangian perturbation expansion in the amplitude of the fluctuations around the 



RW background solution. To this aim, it will prove useful to define the two scalar quantities ( Doroshkevich 1970 ) 



Hi = </? , 7 7 = Ai + A 2 + A3 (32) 



= ^(^o, 7 7 (po, S s -<Po, J s <po, S j ) = A1A2 + A1A3 + A 2 A 3 , (33) 
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where A a are the local eigenvalues of the symmetric tensor (fo, a g. Let us then expand the relevant quantities retaining only 
terms up to second order, 



#°/9 = #(1)0 + $(2)0 
A™f3 = &(l)p + ^-(2)0 
ri — rt (2)/3 



(34) 
(35) 
(36) 



where the zeroth order terms of ^g and vanish, while the expansion of H°j starts from second order. 

For the velocity gradient tensor § a B and the peculiar gravitational field tensor A a g the first order terms obviously coincide 
with the initial data, 

\i)g = -$(1)0 = <Po, a p , (37) 
while the second order corrections read (for more details see Appendix A2) 

V( 2 )0 

A {2)0 = j (20/iiv3o, Q ^ -10n2S a a - 13tpa, a j ipo, J g ) + K a p 

Here indices are raised by the Kronecker symbol. The traceless tensors x"/3 an d K°"g represent the contributions coming from 
the magnetic part of the Weyl tensor to i} a g and A"^, respectively, and have zero divergence like the magnetic tensor TL a g 
itself at this order (see Eqs.(A12)-(A13)). 

From Eqs.(37)-(39) one immediately obtains the traces 



- (-12/11 W0, a f3 +6/i2<5 Q 3 + 5(y30, Q 7 fO^g ) + X% 



(38) 
(39) 



J? = -Mi 
and 



a(-M? + ^Ma) 



(Mi + o(//i 



7 



M2j 



(40) 



(41) 



which coincide with those obtained in the Lagrangian second order Newtonian theory (Buchert 1989; Buchert 1992; Moutarde 
et al. 



1991 



Bouchet et al. 1992; Gramann 1993; Lachieze-Rey 199S). 



We can perform a second-order expansion also for the spatial displacement vector £° 

C = f (0° + £(lf + £(2? 

From Eqs.(25) and (37) one immediately gets 

^ a a j- ol 

Ml) = ~ ai P0, 5(0) 



S(2) = ~ aL P0, 0$ (1) + 

where ij) a g is given by 

V /9 



(2)0 £ (0 ) = a </> /3 5. 



(0) 



3 1 

= - (^^tpo^g +fi 2 S a B + 2<p , a 1 (po^p) + — / dax a 0, 



: £ a and dg" = £(q% one can write 



with trace ijj a a = — |/Lt2- Since cte° 
<Ar a = dg a — aipo, a dq 13 + a'ip a gdq ti . 
Then, the second order perturbative solution for the deformation tensor reads 



T> g = -wpo, 0+a yj 



0> 



(42) 

(43) 
(44) 

(45) 

(46) 
(47) 



we can immediately recognize that the first order term is nothing but the kinematical Zel'dovich approximation (Zel'dovich 
197C0while the second order correction is provided by the symmetric tensor ip a B . 

In Newtonian theory one would write the same formal expression, but the irrotationality condition would automatically 
lead to ip a g — ip, a g, with the potential i[> satisfying the second order Poisson equation 

VV = -|/tt, (48) 

which is consistent with the trace of the relativistic equation (45). This implies that the Newtonian eigenvalues v a of %[)% only 
need to satisfy the condition Va ~ ~ f Z 12 ' wnue m order to get the complete information on the single v^s, one needs the 
Newtonian definition of tj) a g as ip, a g, i.e. a non-local information. The relativistic Va's also solve the Newtonian equations, but 
the reverse is not necessarily true: it just depends upon the specific boundary conditions used in solving Poisson's equation. 
Finally, from Eqs.(26), (29) we can compute the conformal spatial metric tensor 
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a 2 f a 
h a p = S al 3 - 2aifi , al 3 +— (19^0,07 ^o, 7 ^ -12/ii^ ,a^ +6^25 a p) + / daXa/3 ■ (49) 

The traceless, divergenceless symmetric tensor x% representing the contribution to the velocity gradient tensor caused 
by the magnetic part of the Weyl tensor, can be written as a convolution 

X a p(q,a) = J d 3 q'S a e(q')f(\q-ci'\,a) (50) 
of the source tensor 

S a = ^, a p -V 2 (2^ , a p -2^ <p ,\ -<5> 2 ), (51) 
with the time-dependent function /, whose Fourier transform f(k, a) satisfies the third order linear differential equation 

— /+ - — f+ (— + l) —f+ -f - Cx (52) 
dx 3 x dx 2 \x 2 J dx x 

where x = kr, with r = (3/ J 4)t 1/ ' 3 the conformal time, and C = 10A 3 /21fc 4 . The initial conditions are f(xo) = ~£cf( x o) = 
0. 

3.1 Fate of cosmological perturbations 

Let us now analyse the evolution of the solution of Eq.(52), describing the behaviour of perturbations on scales much larger 
(x S> 1) and much smaller (x <C 1) than the horizon scale. 

3.1.1 Outside the horizon: the silent universe 
In the limit x <C 1 Eq.(51) reduces to 

d^ f+ xIx^ f+ x^dx- f ~ Cx ' (53) 
setting g — f + -^x 4 Eq.(53) reads 

— g+-— g + — — 3~0 (54) 
dx :i x dx 2 x 2 dx 

whose formal solutions are g = gox n with n = — 1, —5. Keeping only the growing mode, we then find 

/ „ SL X * = -L A 3 t 4 . (55) 
J 180 378 v ; 

When i« 1, i.e. when the characteristic scale of the perturbation exceeds the horizon distance, x Q /3 ~ (3i 2 /14o)<S < ^, and 
the contribution to $ a p due to the magnetic tensor becomes negligible; a similar reasoning would apply to A^. The relevant 
expressions can then be obtained from Eqs.(38), (39) and (50): 

n « -<Po, a g +^ {-U/jlkpo'p +6^S a fJ + 5<p ,% ^ ), (56) 
A°js ~ <po, a p +^(20^i</?o, a /3 -10//2<5°a - 13(fi , a i Vo, 7 ^), (57) 

njot nj (a 5 T i $ a T "\ v /tro\ 

rt j\ e 7*V°> " + £ f3-yS<fio, <po,v), ■ (58) 

Perturbations with size larger than the Hubble radius evolve as a separate silent universe, with spatial gradients playing no 
role. Note that the expression in Eq.(55), and the related results in Eqs.(56) - (58) could have been obtained by taking the 
c -s- limit in Eq.(52). 

Unfortunately, these general relativistic effects have little implications for the problem of cosmological structure formation, 
since perturbations on super-horizon scales are usually assumed to have very small amplitude, so that a linear approximation 
is sufficient. Conversely, the validity of the silent universe approximation on ultra-large scales might suggest its use as a 
mathematical tool to look for cosmologically interesting background solutions of Einstein's field equations. Nevertheless, there 
are a number of formal consequences of our relativistic solutions, which is worth mentioning. One of these is the absence 
of two-dimensional configurations. This can be seen as follows: if one eigenvalue of <po, a g, e.g. A3, vanishes everywhere, the 
Newtonian theory, with suitable boundary conditions, would imply §3(0) = or xz(a) = qz, i.e. no motion along the third 
axis. This is usually referred as "two-dimensional" gravitational clustering. As far as the second order deformation tensor is 
concerned, one would have v\ + V2 — — § fi2, with /L12 = A1A2, and 1/3 = 0. In the general relativistic case with H a p, instead, we 
immediately find v\ = V2 = —va = — f A*2, and 193(a) ^ for a / ao. The motion dynamically impressed along the third axis 
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soon becomes of the same order of magnitude as that in the other two directions. This effect is to be ascribed to the tide-shear 
coupling term <5°g(i9A — A 7 4 t? 4 7 ) in the evolution equation for the tide, which reduces to — 2/i2<5°g to lowest perturbative order, 
namely second order. The only case when this coupling disappears is when two A a 's simultaneously vanish, i.e. for planar 
symmetry. Therefore $3 (a) = is not an exact solution of the relativistic equations, unless another # a also vanishes. As an 
example, no axisymmetric configurations without motion along the symmetry axis are allowed. 



3.1.2 Inside the horizon: validity 0} the Newtonian approximation 
In the limit x S> 1 it is convenient to define the function g as follows 

" = —/+-/ 

dx x 

where g now obeys the equation 

s_ 

dx 2 



29 + 9 w Cx, 



(59) 



(60) 



with initial conditions at x* ~ 1 equal to g(x+ 
equation reads: 



A 3 xl/54:k 4 and dg/dx(x+) « A 3 xl/18k i . The general solution of the above 



g( x ) = 9~l( x ) + ffa(x) + C [ gi I —xdx-gi 



52 



51 



— xdx 

w 



(61) 



,2 

where g\ and §2 are independent solutions of the homogeneous wave equation -j^ig + g 
Wronskian. Since 

gi — bi cosx, §2 — b2 sinx, W = —6162, 
with 61, 62, integration constants, the general solution reads 
g — bi cosx + 62 sin a: + C (x — x* cos(x — x*) — sin(x — a;*)) . 
Keeping only the growing mode, we find 

Cx 2 2A 3 t 2 



and W 



7-?+-*' 

ax x 



Cx 



f 



21k 2 



If kr 3> 1 we find 
and 



« 5 ~ -3X /3- 

The second order deformation tensor reduces to ip°p = ip, a p, while the conformal spatial metric tensor reads 

h a p = S a j3 — 2a(po, cc j3 +a 2 i(j, CC f3 . 

The remaining relevant quantities reduce to 

v p ~ -fo, p -a<po, 7 <po, fs +2o^, p , 



A a _ a . ( ex. -y 10 / ck \ 

/3 « <Po, p +a I 950, 7 950, '0 -yV) /J ) 



-51* is the 



(62) 
(63) 

(64) 

(65) 

(66) 

(67) 

(68) 
(69) 



All these expressions coincide with those obtained at second order in the Newtonian approximation in Lagrangian form 



(Buchert 1989 



Buchert 1992; Moutarde et al. 1991 



Bouchet et al. 1992; Gramann 1995 



Lachieze-Rey 1993) and can be 



obtained from the c — > 00 limit of Eq.(52). The scalar tp carries information on the influence of the surrounding matter on 
the non-linear dynamics of fluid elements. Note that ip, a p produces a tilt of the principal axes of the first-order deformation 
tensor, <po, a a. 

For the magnetic tensor we find 



H a p + -H° 



0. 



Then, within the horizon TL a p decays as 

H a (a) « ^Ln a s {a„), 
a 



(70) 



(71) 
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Figure 1. Behaviour of the function / vs. x = kr; the dashed lines refer to the two asymptotic expressions in Eqs.(55) and (64), 
respectively for i<1 and x 3> 1. 

where the subscript H refers to the moment when the considered wavelength crosses the Hubble radius. 

A numerical integration of Eq.(52) shows the range of validity of the approximate solutions in Eqs.(55) and (64); this is 
shown in Figure 1, where we plot the numerical solution of Eq.(52) vs. the two asymptotic expressions. 

3.2 Solutions with zero magnetic tensor 

Let us consider now the evolution of an infinite homogeneous ellipsoid, for which TL a a (this can be seen by noting that the 
source tensor in Eq.(51) has no Fourier modes inside the Hubble radius). For this discussion it is more convenient to consider 
as dynamical variables #, A, the shear tensor s a g = ■d"^ — §5"a$ and the tidal one e a g = — |<5°gA. For vanishing magnetic 
tensor, the shear, the tidal field and the metric can be simultaneously diagonalized. The only non-vanishing components of 
s a p and e a p are the eigenvalues s a and e a . In terms of the eigenvalues X a of the tensor ip(, a p we have 

1 i 12 , 5, 2 1 2 10 v 

s a = -A a + -m + a[ — Y^ lXa + 7 A « + 3^1 + 2l^ 2 ) ( 72 ) 

\ 1 j. ( 20 \ 13.2 1 2 26 s 

while the expressions for $ and A are given by Eqs.(40) and (41). 

An alternative expression for the density can be obtained by replacing the volume expansion scalar in the exact continuity 
equation. This will be called the "continuity density" and denoted by a subscript c. We easily find 

3 

l + aA c = (1 + oo/xi) J|(l-aA„ +a 2 v a )~ 1 , (74) 

a=l 

with 

3 6 

Va = ~/i2 — -Aa,(A a -l + Aa+l) (75) 

the eigenvalues of i/) a p (the notation being such that if a = 1, a — 1 = 3, while, if a = 3, a + 1 = 1). Note that A c only 
coincides with A up to second order, i.e. A c = A + 0(a 2 ). 

Let us now give the second order solution for some configurations with degenerate shear and tide eigenvalues, i.e. si = s 2 
and ei = ei- For simplicity, we set ao = in what follows. 

Sphere - In this case A = Ai = A2 = A3, and the eigenvalues of the shear and tide vanish identically, while 

■d = -3A(l + y«A) (76) 

and 
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A = 3A(l + yaA), (77) 

l + aA c = (l-aX- ya 2 A 2 ) 3 . (78) 
Pancake - In this case Ai = A2 = 0, then si = S2 and e\ — e-z- The second order solution is 

si = -ei = y(l+oA 3 ), (79) 
= -A = -A 3 (l + aA 3 ), (80) 
while 

l + aA c = (1-oAs)" 1 , (81) 
coincides with the exact comoving density (cfr. Matarrese, Pantano & Saez 1993). 

Filament - In this case we have Ai = A2 = A and A3 = 0. Again, the shear and the tidal tensors have two equal eigenvalues 
which at second order read 

ai = «2 = ~(l-yoA), (82) 

Cl=C 2 = |(l-yOA), (83) 

while 

■& = -2A(l - yaA) (84) 
and 

A = 2A(l + yaA), (85) 
l + aA c = (l-aA- y« 2 A 2 ) 2 (l + ^a 2 \ 2 ^j \ (86) 



4 DISCUSSION 

Our second order perturbative approach to the general relativistic dynamics of irrotational dust can be used to provide some 
insight on the main dynamical issue: the fate of general perturbations. So far, only a few analytical solutions of the Lagrangian 
system of Eqs.(15) - (24) are known: for exact planar symmetry, Ai = A2 = 0, one recovers the Zel'dovich pancake solution 
(Zel'dovich 1970), as shown by Matarrese, Pantano & Saez (1993), which, in a relativistic context can be seen as a particular 
Szekeres model (Szekeres 1975; see also the discussion by Kasai 1992, 1993 and Croudace et al. 1994); in the spherically 
symmetric case, Ai = A2 = A3, the system admits the so-called Tolman-Bondi solution (Tolman 1934, Bondi 1947; see also 
the discussion by Matarrese, Pantano & Saez 1993), which, from a local point of view, has a direct Newtonian analog in the 
well-known top-hat model (e.g. Peebles 1980). A systematic study of the local system of equations (i.e. that obtained by 
setting the magnetic tensor to zero) has been carried out by Croudace et al. (1994), who looked for solutions representing 
attractors among general trajectories. They found that both spherical collapse and a perfect pancake are repellers for general 
initial conditions, but argued that the pancake instability would probably disappear in the presence of a non-zero Ti a p- 

Bertschinger & Jain (1994) have shown that the instability of the pancake solution in the local system of equations is 
caused by the tide-shear coupling in the evolution of the tide, which has a destabilizing effect on the pancake solution (for 
general initial conditions) but stabilizes prolate configurations. According to this analysis, for vanishing H a p, a strongly prolate 
spindle with expansion along its axis is the generic outcome of collapse, except for specific initial conditions corresponding to 
exactly spherical or planar configurations. They proposed a Newtonian interpretation for this collapse theorem and argued that, 
even though in the most general case 7i a p is likely to be non-zero, it can be set to zero for highly symmetrical configurations 
and probably neglected in many other circumstances. 

As stressed by Matarrese, Pantano & Saez (1994), the tide-shear coupling term in the tide evolution equation is not 
present in the equations of Newtonian theory, at least in its standard formulation, although it is clearly compatible with 
them. The dominant role of this term for the non-linear fluid evolution is a peculiarity of the purely local equations, which 
certainly applies to perturbations on super-horizon scales, but cannot be easily understood in the Newtonian limit. As a 
further illustration of this interpretation let us consider the collapse of an infinite homogeneous ellipsoid, which is described 
by our equations when H a p = 0. In such a case the Newtonian dynamics is known to lead to the formation of oblate spheroids 
(e.g. White & Silk 1979; Barrow & Silk 1981), pancake-like objects with one collapsing axis and the other two tending to 
a finite size (apart from initial conditions corresponding to an initially prolate spheroid). The general relativistic dynamics, 
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instead, favours the formation of prolate objects, collapsing filaments with expansion along their axis. This point is discussed 
in some detail by ZePdovich & Novikov (1983), who argue that the discrepant behaviour is due to the different role of 
boundary conditions at infinity in the two theories. Our second order calculations show that the evolution of fluid elements 
as isolated ellipsoids cannot apply to perturbations on scales smaller than the Hubble radius, where non-local effects would 
play a fundamental role. The appearance of non-local terms in general solutions of the Newtonian theory has been recently 
stressed by Kofman (1994). 

Given that the magnetic part of the Weyl tensor is generated during the mildly non-linear evolution of irrotational dust, 
even in case it was primordially set to zero, is there any application left for a purely local treatment on cosmologically relevant 
scales? One can try to make some guess. It may be that during the late phases of the evolution of a given fluid element, 
i.e. in the strongly non-linear regime, the role of the magnetic tensor becomes negligible again, possibly due to the fact that 
some sort of (approximate) local axisymmetry is dynamically established. However, if this picture should prove correct one 
cannot assume that the principal axes of the collapsing ellipsoid coincide with the initial ones; in fact, already at second 
order in perturbation theory the effect of the magnetic tensor is that of tilting these axes, on account of the interaction with 
the environment. A further complication is induced by the occurrence of orbit crossing with the subsequent formation of 
multi-streaming. 

Let us finally mention what we think is an important point. The main lesson we have learned from the purely local 
treatment, our silent universe approach, is that one could try to isolate two competing effects: i) If H a p is switched off, the 
fate of the collapse of a general fluid element is completely determined by its local initial conditions; in such a case, one would 
generally expect prolate rather than oblate collapse, ii) In the real world, however, H a /s is non-zero except for a number of 
unrealistic cases listed in the introduction; H"^ carries information on the action of the "rest of the world" on the considered 
fluid element, i.e. what can be called the effects of the environment. In other words, the matter surrounding the fluid element 
acts by compressing or stretching, and generally deforming it, thus modifying its local dynamics in a way which cannot be 
easily predicted a priori. The relevant issue is: which of these two competing effects is going to dominate on the long run? It 
seems to us that the answer generally depends upon various variables: the scenario (e.g. cold or hot dark matter, amount of 
baryons, etc.), and the overall initial conditions, as they are usually implied by the power-spectrum and the statistics of the 
primordial density fluctuation field. 

While completing this work two preprints have circulated (Bertschinger & Hamilton 1994; Kofman & Pogosyan 1994), 
where the non-linear evolution equation for the electric component of the Weyl tensor is derived from Newtonian gravity. The 
presence of a non-vanishing non-local contribution in this limit agrees with our previous results (Matarrese, Pantano & Saez 
1994) on the behaviour of the general relativistic equations inside the Hubble radius. 
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APPENDIX A: NEWTONIAN DYNAMICS 

The equations which govern the non-linear dynamics of a collisionless fluid in Newtonian theory for an expanding universe 
can be written as (e.g. Peebles 1980) 

^ + (l+<5)V-v = 0, (Al) 

5 + 2-v+4rV0 = O, (A2) 
dt a a 2 

V 2 <j> = ^a 2 8. (A3) 

where v = dx/dt is the peculiar velocity, (f> the peculiar gravitational potential and 5 = (6irGt 2 p — 1) the fractional density 
contrast. In order to compare these equations with the relativistic ones it is convenient to use rescaled variables for the velocity 
field u = da/da = v/d, the density fluctuation field A = 5/a and the local gravitational potential if = (3/2A 3 )(f> and to adopt 
as time variable the scale factor a itself (e.g. Gurbatov, Saichev & Shandarin 1989; Matarrese et al. 1992). Eqs.(Al-A3) in 
the new variables read 

A' + ^A, /3 =-i(«' 3 ,/3+A)-A u %, (A4) 

u°" + u?u a ,p=-?-(u a + V , a ), (A5) 

<p,%=A. (A6) 

By differentiating the Euler equation (A5), defining the symmetric tensors ^"g = u a ,p, A a /3 = (p, a p, and adopting a Lagrangian 
description, one recovers Eq.(15), while the continuity equation (A4) coincides with the trace of Eq.(16). It is clear that the 
Newtonian approximation described above is degenerate, as it provides only one equation to fully determine the symmetric 
tensor A a g : any traceless tensor added to the r.h.s. of Eq.(16), leaves the Newtonian equations unchanged. In order to 
completely determine the evolution of the gravitational field tensor one has to resort to its definition in terms of the 
potential <p, i.e. to a non-local theory. Because of the intrinsic non-locality of the Newtonian approximation one needs 
boundary conditions to determine the dynamics, contrary to the general relativistic equations, where Cauchy data on a 
spatial hypersurface are enough; this comes from the fact that the Poisson equation (A6) is an elliptic, constraint equation. 
It is well known (e.g. Ellis 1971) that the lack of evolution equations for the traceless part of the gravitational tensor, e"^, 
implies that the Newtonian theory may add spurious solutions which would be discarded by the full relativistic system. 



APPENDIX B: SECOND ORDER SOLUTIONS 

At second order the differential equations (15)-(17) become (for convenience the integration variable is now the proper time 
of fluid elements) 

#W + 7? faw + A (2 fc) + df , a -, W?p = 0, (Bl) 



2a 

A(2)g + ^ ($(2)0 + A (2 )g) + a (3v?o, a 7 <fio, J -4/xnpo,°0 +2/i 2 <5 Q 3 ) = ^- (eJ s H a 7 ,s +e c " ri 6 H"' ,5 ) , (B2) 



' a 



2a \ 



e« A (2)7 ,(S +e -v A 



7 

(2)01 



5A 3 , 



7^ 



<P0, 



cvyS 



) <fio, 



(B3) 
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where the last term in Eq.(A9) has been presented for completeness, but it is a decaying mode and will be neglected in what 
follows. Note that in order to write Eqs.(A8) and (A9) we used the Christoffel symbols of the three-metric h a p evaluated at 
first order; these are trivially obtainable in terms of the initial gravitational potential tpo and its derivatives. 

In the expressions for $(2)p and A(2)p i s convenient to isolate those parts (indicated by a tilde) which solve the system 
(A7)-(A9) with vanishing magnetic term H a p. It is possible to proceed in this way because the contributions from the magnetic 
term start only from second order. We have 

*we = *w + x% (B4) 

A (2 )0 = &(2)g + K a g. (B5) 

The traceless tensors \ a p and represent the contributions caused by the presence of the magnetic part of the Weyl tensor 
and have zero divergence, 

Xlt;a = Kp" :a = 0. (B6) 

Similarly, at second order, we have 

H(? ia = 0. (B7) 

The previous conditions are directly derived from the constraint equations, and are obviously consistent with the evolution 
equations. 

For zero magnetic term, the solutions of Eqs.(A7) and (A8) can be easily found; they read 
= |(-12Mi^o, a /3 +eH2d a g + 5<p ,% <po,\ ), (B8) 

A (a fr = j(W^ lVo , a p -W^S a p - %,>o,;). (B9) 
The second order terms due to the magnetic tensor satisfy the following equations 

X a + Y a + ^) = °' ( B1 °) 

* a + I (X a + K a p ) - Ya ( e / W V +£ a74 Wfl7.« ) = 0, (Bll) 

H a p + -H a p = ~ (c/kV +e a T Kfh , f ) - y {C a p + Cp a ) , (B12) 
where 

C a g = e" 7 * (<po, p y <po,Sv ) " ■ (B13) 
One can easily eliminate 7i a p by differentiating Eq.(A17) with respect to t. One gets 

-a - 2 . a 2 . a 2/ a a \ ^ To ck7?7 u<5 . uS ( 'yri a . a"/ri \ 1 

K + J n + g^X + 7^2 [X + K 0) = [ 2e £ K ~'HM + e 7 (V K + £ K0H,S V )\ 

~ [e/ {C\ s + C^ s ) + e^' 6 (C 7/M + Cp^)} . (B14) 



Using the relation 

e a/37 W = 3! StyW, (B15) 
we then obtain 

« + -« - [X + k 0) - k = -— S p, (B16) 
where the constant source tensor S a p is given by 

S a p = V2, a p -V 2 (2^ , a p -2^ , a 7 ^p -5> 2 ). (B17) 
Replacing now 

* a = -| (^) (B18) 
into the differential equation (A22) we get 

d 3 a 5 d 2 a 32 Ci a 1 2 / ^ a 1 a \ 10 Q /Dln > 

d^+id^+K** 3 -^ \dt x "--t x n = -^-t s ^ (B19) 

which has to be solved together with the initial conditions 
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x a (to) = j t x%{t ) = ^x a (to) = 0, (B20) 

since our linear initial data imply X a /3( a o) = K °/3( a o) = H. a p(cio) = 0. 

The traceless tensor x"p, representing the contribution due to the magnetic part of the Weyl tensor, can be written as a 
convolution 

X Q /3 (q,a) = ydV^(q')/(|q-q'|,a) (B21) 

of the source S a p with the function /, whose Fourier transform f(k,t) satisfies the equation 
d 3 - 5 d 2 f 32 d f k 2 ( d f 1 A _ 10 

dtS f+ tdtt f+ WTt f+ ^{Tt f + l f ) - Tat' (B22) 

with initial conditions f{to) = [d//dt] (to) = [d 2 //dt 2 ] (to) = 0. Note that the possibility of 3D Fourier expanding the 
previous quantities is ensured by the underlying flat RW background. 

In order to discuss the general behaviour of the previous equation it can be useful to adopt as independent variable 
x = kr where r = (3/A)t 1//3 is the conformal time. One then gets Eq.(52) whose properties are discussed in the main text. 
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